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Abstract 



We use a multitype continuous time Markov branching process model to de- 
p { scribe the dynamics of the spread of parasites of two types that can mutate into 

each other in a common host population. Instead of using a single virulence char- 
acteristic which is typical of most mathematical models for infectious diseases, 
our model uses a combination of two characteristics: lethality and transmissi- 
bility. This makes the model capable of reproducing the empirically observed 
fact that the increase in the host density can lead to the prevalence of the more 
virulent pathogen type. We provide some numerical illustrations and discuss the 
effects of the size of the enclosure containing the host population on the encounter 
rate in our model that plays the key role in determining what pathogen type will 
eventually prevail. We also present a multistage extension of the model to situa- 
\ tions where there are several populations and parasites can be transmitted from 

' one of them to another. 



1 Introduction 

There is an extensive literature on mathematical modeling of infectious diseases. Both 
the rate of spread and the virulence of pathogens are important, and both the dynamics 
of disease spread and the evolution of virulence in parasite-host systems have attracted 
much attention. The first mathematical model of disease dynamics apparently goes 
back to the 1760 D. Bernoulli's paper [13] (see also pp. 2-6 in [16]). A major step 
in further development of deterministic modelling was the "threshold theorem" [38] 
showing that the initial frequency of susceptibles must exceed a certain threshold value 
to give rise to an epidemic. The first stochastic epidemiological models appeared in 
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the late 1920s [43], applying the "law of mass action" suggested in [13] to probabilistic 
description of the epidemics. Stochastic models became much more popular after 
Bartlett [8] formulated a model for the general stochastic epidemic by analogy with 
the deterministic model from [38]. May and Anderson [47] introduced the first model 
of evolutionary change in pathogen virulence in 1983, which explained the decline 
in virulence of myxomatosis, introduced to control Australian rabbits, in terms of 
optimal transmission of the parasite between hosts. Subsequent work, notably [47, 
19, 44, 40, 23], has expanded on this classic paper, to cover the effects of a vector 
of the pathogen, vaccination of hosts, and various other issues. For more extensive 
surveys of the field and literature reviews, we refer the reader to the monographic 
literature [2, 3, 6, 10, 16, 49] and also to papers [4, 18, 34]. 

In natural systems, the density of host organisms declines when a virulent pathogen 
produces an epidemic that kills host rapidly. The models of virulence have made clear 
that this leads to selection for less virulent strains, such that infected hosts live longer, 
and thus tend to pass the pathogen to more susceptible new hosts, maintaining the 
epidemic. While hosts will have longer generation times than their pathogens, they 
too, will evolve to be more resistant, further reducing the virulence of the pathogen for 
these hosts. Thus one can expect pathogens to be highly virulent only when they have 
recently switched hosts, such as the SARS coronavirus that was probably transferred to 
humans from animals such as Himalayan civet cats [29] or when the plague bacterium, 
Yersina pestis, has been transferred to humans from rodents [56]. 

The main objective of the present paper is to provide quantitative and also quali- 
tative (cf. the discussion of the above-mentioned threshold theorem and other findings 
from [38] on pp. 11 and 29 in [16]) insights into the dynamics of pathogen populations 
where the virulence of the pathogens can change due to mutation. More specifically, 
we are interested in analyzing the effects that changes in host density may have on 
the virulence of the pathogen. In agricultural systems, humans have been maintaining 
animals and plants at very high densities for about 10,000 years, and humans them- 
selves began living at high local densities in villages, towns and cities from about the 
same time [15]. Densities of both humans and domesticated organisms have probably 
increased continually since then, but the densities at which some domestic animals 
are kept appear to have increased very rapidly in the last few decades (see e.g. [22], 
and the recent rapid development of large scale commercial aquaculture involves the 
maintenance of very high densities of a whole new set of marine species that will have 
their own pathogens (see [32, 31, 12]). Thus the question arises as to whether high 
host densities are likely to lead to selection for more virulent strains of pathogens. 

The mathematical models we discuss in the present paper incorporate a continuous 
time Markovian multitype branching process. Branching processes first appeared in 
epidemiological context in the mid-1950s: one can mention here the paper [14] and 
also refer to [9, 37, 55]. For further literature review, we again refer the reader to the 
above-listed monographic and survey literature, while some examples of applications 
of multitype branching processes in epidemiology can be found in [11, 48, 33]. 

Branching processes are relatively simple and well-understood mathematical mod- 
els; for treatments of the theory of branching processes, see e.g. [35, 5], a recent addition 
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to the monographic literature on branching processes prepared specifically for biolog- 
ical audience being [30]. At the same time, the processes were shown to be good 
approximations to the general stochastic epidemic models at the initial stage (when 
the total population size is large and the initial number of infectives is small) and also 
at the final stages of the epidemics (see e.g. [42, 7, 34]). As our main objective is not 
to model the detailed dynamics and describe the whole history of the epidemic itself, 
but rather to suggest a model for changes in the pathogen's population composition 
depending on the density of the host population, we will restrict ourselves to working 
with the simpler branching process models. 

Moreover, the most critical stage of an epidemic is the initial one, when it is basically 
determined if there will be a large-scale event or the epidemic will die out. And as 
branching processes are good approximations to the general stochastic epidemic models 
at the initial stage, the threshold analysis aimed to determine if the "basic reproductive 
number" (defined roughly as the expected number of secondary cases produced by one 
infected individual) is greater than one (which implies the danger of an outburst or 
persistence of endemic levels) can be carried out using those models (see e.g. Chapters 6 
and 8 in [49]). 

An additional argument for this approach is that in agriculture or aquaculture 
situations, when a threat of an epidemic outburst arises, drastic measures are usually 
taken: the affected part of the population may be isolated or even destroyed. These 
measures will substantially change the way pathogens can pass between hosts, and 
hence can make the standard stochastic epidemic models inapplicable beyond the initial 
stage of the epidemic anyway. 

One can also envisage an extension of our simple analysis to the general stochastic 
models as the same mechanism will certainly work for the latter as well. However, 
such extensions will be much less tractable analytically and may lead to no closed-form 
answers. 

In our analysis, we will look at supercritical two-type branching processes (so that 
the basic reproductive number will be greater than one: we are interested in what 
happens when there is a threat of epidemics) and then look at the behaviour of the 
ratio of the sizes of the subpopulations in the process (representing two versions of the 
pathogen, that can mutate into each other). This quantity can be used to determine 
which of the two types will become dominant in the population over time. One of the 
advantages of the branching process model is that one can easily incorporate in it two 
different characteristics of pathogenes' virulence: their lethality (which can be described 
by the mean survival time of an infected host) and transmissibility (specifying the 
probability of an infected host infecting a susceptible one on their contact), instead of 
representing them together by a single quantity termed simply "virulence" (cf. e.g. [1]). 
For recent discussions of the interplay between these two characteristics, see [27, 41] . 

In the paper, we use our approach to model the dynamics of a host /parasite popu- 
lation where parasites can be of one of two types that can differ in their lethality and 
transmissibility. The underlying simple continuous time Markov model of a two-type 
branching process is presented in Section 2. The analysis of the model and derivation 
of the dynamics for the mean functions are given in Section 3 and used in Section 4 
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to establish the eventual composition of the pathogen population. It turns out that 
our model is capable of reproducing the above-mentioned phenomenon of shifting the 
overall parasite's lethality in response to increased density of the hosts. This is actually 
possible because of using two different parameters to characterise the virulence of the 
pathogen. 

In Section 5 we present a few remarks on how the change in the size of the enclosure 
a given population of hosts inhabits can affect the "encounter rate" for the hosts — the 
key parameter of the model describing the "effective density" in the host population. 
Finally, in Section 6 we consider a multi-stage modification of our model that can be 
used to analyse farm or aquaculture situations in which there are many enclosures or 
tanks of animals, and an outbreak of an infectious disease occurs in one of them. We 
assume that the pathogen might at some stage be transmitted to the next enclosure, 
where the epidemic process starts anew etc., and discuss possible scenarios of the 
development of the epidemic in the farm. Section 7 contains a few final remarks 
concerning biological interpretation of our results. 

2 Description of the branching process model 

Assume that we have a large population of hosts that can be infected by parasites of 
one of two types that will be denoted by T\ and T 2 . The pathogen types can differ in 
both their lethality and transmission rate. The numbers of infected hosts at time t are 
represented by the vector 

Z(t) = (Z 1 (t),Z 2 (t)), t>0, 

Zi(t) being the time t number of hosts infected with the type T« pathogen (for simplicity, 
at any given time, any given host can be infected with one type of the pathogen only). 
We do not keep track of the number of hosts that remain uninfected (susceptible), 
assuming instead that this number will remain large enough during the time period 
for which our mathematical model is intended, so that the dynamics of the process 
{Z(t), t > 0} to be described do not change over time. 

A general description of the model is as follows. We assume that each infected 
individual lives a random time (which will tend to be shorter when one is infected 
with the "more lethal" of the two pathogen types). During its lifetime, an infected 
host can encounter susceptible hosts and, with a probability depending on the type 
of the pathogen it carries, transmit the parasite to them. The rate of such (random) 
encounters will be specified by a special parameter that we can vary in order to model 
changes in the density of the host population. 

Finally, we also allow the pathogens to mutate, so that when a host originally 
infected with 7\ encounters a susceptible host, the latter can become infected with 
2V type parasites (and the other way around). 

Now we will present a formal mathematical model. First recall that the exponential 
distribution with rate (or intensity) a > has density of the form 

p(t) = ae~ at , t > 0, (1) 
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with mean 1/a, and plays a special role in probability theory due to its unique memory- 
less property that makes the distribution ubiquitous in the theory of Markov processes. 
Namely, a random variable r > modelling, say, the time of the first encounter of a 
given infected host with a healthy one, has this property if, for any s,t > 0, the 
conditional probability of the event {r > s + t} given that r > s coincides with the 
probability of {r > t}: 

iw , *i \ Pr(r > s + t, t > s) Pr(r>s + t) 

Pr(r > s + t\r > s) = — - — = — ±-. —- = Pr(r > t). (2) 

V 1 ; Pr(r > s) Pr(r > s) V > \ ) 

In words, if, at time s we know that there has been no such encounter, then the 
conditional distribution (given that information) of the residual random time r — s till 
the encounter will be the same as the original distribution of r. It is obvious that if r 
has density (1) then 

/oo 
p(s)ds = e~ at , t > 0, (3) 

and so the property (2) is clearly satisfied. 

An equivalent formulation of the property can be given in terms of the distribution's 
hazard rate r T (s) that quantifies the probability that, given there has been no encounter 
up to time s, there will be one "immediately afterwards": in case a random variable r 
has a continuous density p T , the hazard rate is defined by 

r ^ = P^T7)--S kPr(T> * s>0 - (4) 

and, as h i 0, 

Pr(r < s + h\r > s) = r T (s)h + o(h), s > 0, (5) 

where, as usual, o(h) denotes a quantity that vanishes faster than h: o{h)/h — > 0. 

It is obvious from (3) and (4) that the hazard rate of a distribution is constant if and 
only if it is exponential (in that case, the hazard rate simply equals the distribution's 
rate). In applications, one often uses exponentially distributed random variables to 
model times between successive events of a particular kind and also lifetimes. This is 
because, on the one hand, such assumptions make sense from the modelling view point 
(in a large population, meeting a new individual during a time interval [t, t + h], h > 0, 
can scarcely depend on one's "life history" prior to time t) and, on the other hand, 
as the resulting models are usually Markovian, so that the powerful machinery of the 
theory of Markov processes is applicable. 

Our basic model assumptions are as follows: 

(a) The initial population contains Zi hosts infected with parasites of type Tj, i — 1,2. 

(b) A host infected with type Tj parasites lives a random time which is exponentially 
distributed with parameter > 0, i = 1,2. Clearly, the pathogen with a higher 
rate will be the more lethal one, as the mean lifetime in that case will be lower. 
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(c) Any infected host can encounter susceptible ones. The time till the first encounter 
of a given infected host (of any type) with a susceptible host is a random variable 
exponentially distributed with rate A. It is clear from the above discussion of the 
properties of exponential distributions that, at any time t, the residual times till 
encounters of the current infected hosts with susceptible hosts are all exponential 
with the same rate A. A similar observation applies to the residual lifetimes; this 
ensures that the process {Z(t)} will be Markovian: given the current state of the 
process, its future evolution does not depend on the past one. 

In a modification of the model, one can assume that, at time t, the time till the 
first encounter of a given infected host (of any type) with a susceptible host has 
a hazard rate X(t) which can depend on t. This will enable one to model changes 
in the density of the host population that occur over time (the higher A, the more 
often are encounters, which corresponds to higher host density situations). The 
process will remain Markovian, but will become time-inhomogeneous. 

(d) At any encounter with susceptible hosts, a Tj-infected host meets only one sus- 
ceptible host (there can be several such encounters during the host's lifetime). At 
each such instance, the T r infected host transmits the parasites to the susceptible 
one with probability fa (so that, with the complement probability 1 — fa, the 
encounter will have no consequences for the susceptible host). 

(e) Mutations Ti -h- T 2 are possible. A host infected with Ti-type pathogens will 
remain such till its end, but, when transmitting pathogens to a susceptible host 
during an encounter, the newly infected host will carry T 2 -type pathogens with a 
probability H\. Likewise, /z 2 denotes the probability that a successful transmission 
of parasites from a T 2 -infected host to a susceptible one resulted in making the 
latter Ti-infected. 

(f) All the above-mentioned random times (lifetimes, times till the first encounter) 
are independent of each other. 

Of course, the above assumptions oversimplify the real biological processes. There 
are likely to be several or many strains of pathogen, and the probability of mutation 
from one to another will vary. We suggest that simplifying the system to two strains is 
likely to retain the same key dynamic features. Further, the distributions of times until 
events occur are likely to be approximately exponential as argued earlier, and we do 
not see any reasons why host survival times and encounter rates should depend in any 
way on each other. Note that "encounters" between hosts are simply occasions where 
a pathogen can be transferred between hosts. They do not have to involve physical 
contact. Thus a pathogen transferred by aerosols might be transferred between pigs 
that are isolated in neighbouring stalls. 

The diagram in Fig. 1 illustrates the "physical" meaning of our assumptions. The 
two horizontal segments represent the mean lifetimes of hosts infected by the pathogens 
of our two types: the longer segment (of length l/cti) corresponds to Ti which we 
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Figure 1: Mean lifetimes of infected hosts and mean times to encounter. 

assume less lethal by stipulating that a 2 > oti, whereas the shorter one (of length 
1/0:2) corresponds to T 2 . 

When the host population density is relatively low (say, represented by the value 
A = A", as depicted), the lifetime of T 2 -infected hosts will be too short compared to 
the mean time 1 / A" between encounters to give them a good opportunity to encounter 
susceptibles and hence further propagate in the host population. One may expect that 
this will result in the eventual prevalence of 7\ pathogens who have better chance 
of being transmitted as they live longer. However, if the host population density 
increases (say, to A = A', as depicted), then the more lethal type T 2 may have frequent 
enough encounters which, combined with its higher transmissibility, can lead to its 
eventual prevalence. As we will see in Section 3, the above argument is confirmed by 
mathematical analysis. 

Assumptions (a)-(f) imply that our process {Z(t)} is actually a two-type time 
homogeneous Markov branching process in continuous time, see e.g. Section V.7 in [5]. 
That is, {Z(t)} describes the dynamics of a population consisting of individuals (or 
"particles") of two types, 1 and 2. The transitions of different particles in the process 
are assumed to be independent. A particle of type % lives for an exponentially random 
time with rate a,. At the end of its life it disappears. It can either 

(i) simply disappear (in terms of our modelling assumptions above, this means that 
a given Tj infected host died having never encountered a susceptible host), or 

(ii) produce one particle of the same type (meaning: there was an encounter, but 
no transmission occurred; we think of the "newly produced" particle of type i 
as just the "old" infected host of type T who keeps living — note that, due to 
the memoryless property of the exponential distribution, such an identification 
of a "new" particle with the "old" one is in agreement with our assumptions 
(a)-(f) (so that the life of one Tj-infected host can actually be represented by a 
succession of several type i particles, for which reason we do not use Tj to denote 
the type of particles in the branching process), or 

(iii) produce two particles of the same type % (meaning: there was an encounter and 
successful transmission, but no mutation; one of the "newly produced" particles 
is actually the original host, the other represents a newly infected — with the 
same Tj-type pathogen — host), or 
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(iv) produce two particles of different types (meaning: there was an encounter and 
successful transmission and mutation; one of the "newly produced" particles 
is actually the original host, the other represents a newly infected — with the 
pathogen of the other type — host). 

To see that both sets of assumptions (a)-(f ) and (i)-(iv) describe the same dynamics 
of {Z(t)}, it suffices to note that in both cases we deal with time- homogeneous contin- 
uous time jump Markov processes (which follows from the exponent iality assumptions) 
and then verify that, choosing suitable parameter values for the second model, one can 
obtain the same transition rates as for the first one. 

To do that, we first observe that assumptions (i)-(iv) imply the branching property 
which means that, for any s > 0, given Z(s) = (z±, z 2 ), the future {Z(s + £),£> 0} of 
the process will follow the same probability laws as that of the sum of z\ independent 
copies of Z(t) starting at time with a single particle of type 1 and z 2 independent 
copies of Z(t) starting at time with a single particle of type 2. It is clear that the 
first model (specified by (a)-(f)) has the same property. 

Moreover, the branching property implies that to completely describe the evolution 
of the process, it suffices to specify transition probabilities 

P?\jnh) = H Z ( h ) = Ui,h)\ Z (°) = e 
from the basic initial states 

ei = (l,0) and e 2 = (0, 1) 

for arbitrary small time increments h. It is obvious that the probability of having more 
than one transition during a small time interval (0, h) will be o(h), so we just need to 
consider where a single transition can take the process from a basic state e$ according 
to assumptions (a)-(f) and show that the transitions will have the same rates as for a 
process specified by (i)-(iv) (for a suitable choice of parameter values). 

Suppose that, in our branching process, particles of type % have exponentially dis- 
tributed random lifetimes with rates = on + A, % = 1,2. Moreover, at the end of a 
type % particle's life, it produces a random number of children (possibly of both types) 
according to the offspring distributions qi(ji,j 2 ) = Pr(a particle of type % gives birth 
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to ji particles of type 1 and j 2 ones of type 2) given by the following table: 



U1J2) 


Q1U1J2) 


Q2U1J2) 


(0,0) 


a\ 

A + «! 


a 2 
X + a 2 




(l-A)A 
A + «i 







(l-//i)/3iA 
A + «i 


n 
u 


(1,1) 


/iiftA 
A + ai 


/^AjA 

A + a 2 


(0,2) 





(1 - ^ 2 )/5 2 A 
A + a 2 


(0,1) 





(1-/3 2 )A 
A + a 2 



Table 1. Offspring distributions in the branching process model (i)-(iv). 
Then, due to independence and (5), 

p^(0,0) = Pr(initial type 1 particle dies in (0, h), no children produced) 
= Pr (initial type 1 particle dies in (0, h)) x q 1 (0, 0) 

= ((cti + X)h + o(h)) x = aih + o(h), 

A + «i 

which is clearly the same as the probability for a single Ti-infected host (from the first 
model) to die during (0, h). 
Likewise, 

p^(2, 0) = Pr (initial type 1 particle dies in (0, h), producing 2 children of type l) 
= Pr (initial type 1 particle dies in (0, h)) x q\(2, 0) 

= ((«! + X)h + 0(h)) X t 1 -^^ = (1 _ ^)^\h + o(/l). 

A + Qt\ 

The corresponding transition in the first model is as follows: a Ti-infected host did not 
die during (0,h), but met a susceptible; the pathogen was transmitted, no mutation 
occurred. Due to independence, the probability of this will be 

(1 + 0(h)) x (Xh + o(h)) x ft x (1 - /Xi) = (1 - //O&A/i + o(h). 

Virtually the same argument shows that 

Pi^(l, 1) = Pr (initial type 1 particle dies in (0, h), producing one child of each type) 
= Pr(initial type 1 particle dies in (0, h)) x q±(l, 1) 

= ((^ + X)h + o(h)) x = j J , 1 j3 1 Xh + o(h), 

A + «i 
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coincides with the probability for a Ti-infected host to meet a susceptible and transmit 
the mutated form of the pathogen. 

Finally, given Z(0) = e±, the most likely state for Z(h) is e.\. In our branching 
process this occurs when either the original particle lives through the time interval (with 
probability 1 — (a.\ + X)h + o(h)) or it dies prior to h producing a single type 1 particle 
(of which the probability is ((a 1 + \)h + o(h)) x (1 - /3 1 )X/(X + a 1 ) = (l-/3i)Xh + o(h)), 
so that 

p[ h \l, 0) = 1 - (e*i + X)h + (1 - fi^Xh + o(h) = 1 - (an + /3iX)h + o{h). 

In the first model, to get to this state, the initial Ti-infected host survives for h time 
units and either has no encounters with susceptibles or has one, but without transmis- 
sion of a pathogen. The probability of this is 

(1 - ai h + o{h)) x [(1 - Xh + o{h)) + ((1 - fa)Xh + o{h))\ 

= (1 - a x h + o(h)) x (1 - f3 ± Xh + o(h)) = 1 - (a x + p x \)h + o(h), 

the same value as for the branching process model. 

Of course, we could also work out the probability for (1, 0) just as the complement 
probability 

^(1,0) = 1 - (rf fc) (0,0) + P P(2,0) +ri fc) (l,l)) +o(h), 

but the presented argument demonstrates the difference between the interpretation of 
the elements of the two models (Tj-infected hosts in the first model and type % particles 
in the second one are not the same, as we noted earlier). 

The same calculations are applicable in the case of the initial state e 2 , which leads 
us to the transition probabilities presented in Table 2 (for h — > 0, the additive terms 
o(h) being omitted for brevity). 



U1J2) 






(0,0) 


aih 




(1,0) 


1 - (c*i + (3iX)h 





(2,0) 


{i-m)PiXh 





(1,1) 


/iiPiXh 


liifiiXh 


(0,2) 





{l-^ 2 )(5 2 Xh 


(0,1) 





1 - (a 2 + (3 2 X)h 



Table 2. Transition probabilities for small h values (o(h) terms omitted). 

3 The dynamics of the means 

Consider the matrix of mean values M(t) = (M^ (t)), where 

M i5 {t) = Z(0) = ei ), i,j = 1,2, 
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is the expected number of type j particles present in the process at time t given that 
the process started at time with a single particle of type i. Clearly, 



M(0) 



1 
1 



the identity matrix, and, using the branching and Markov properties of {Z(t)}, it is 
easy to demonstrate that {M(t), t > 0} possesses the operator semigroup property: 



M(s + t) = M(s)M(t) for any s, t > 0. 



(6) 



Indeed, given that Z(0) = (z 1 , z 2 ), the value of Z(s) is just the sum of z\ independent 
copies of Z(s) starting at time with a single particle of type 1 and z 2 independent 
copies of Z(s) starting at time with a single particle of type 2. As the process is 
time-homogeneous, we infer that 



E(Z(s + t)\ Z{s) = ( Zl ,z 2 )) =^2ziE(Z(t)\ Z(0) = ei) = ( Zl ,z 2 )M(t). (7) 



i=i 



From here, using the Markov property and the double expectation law for conditional 
expectations, we have 

E(Z(s + 01 Z(0) = ei) = E[E(Z(s + t)\ Z{s)) | Z(0) = e % ] 

= E[Z(s)M(t)\ Z(0) = e t ] = (M^s), M l2 (s))M(t), 

which is equivalent to (6). 

Relation (6) implies (cf. p. 202 in [5]) that one has the matrix exponential repre- 
sentation 

oo ,k Ak 



k=0 



where A = I and 



A = lim-(M(t) - I) 

Ho h 



(9) 



is the so-called infinitesimal generator of {M(t)}. Indeed, from (6) one has M(t + h) — 

M(t) = {M{h) - I)M(t), so that (9) implies that jM(t) = AM(t), M(0) = /, for 

which (8) is clearly a solution, as seen from its term-wise differentiation. 

Evaluating matrix exponentials is rather straightforward: it basically reduces to 
calculating the values of the function on the matrix' spectrum (for more detail on 
functions of matrices, see e.g. Chaper V in [24]). If, say, A is diagonalisable, so that 
there exists an invertible matrix Q (with inverse Q l : Q l Q = QQ 1 = I) such that 



Q AQ = D = diag {a + , a_} = 



a+ 

<7_ 



(10) 
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for some a± (which is the case in our situation, as we will see below), then clearly 
A = QDQ \ 

A 2 = (QDQ l ) 2 = QDQ ^QDQ 1 = QD 2 Q 1 = Q diag {a 2 + , a 2 _} Q 1 
and so on: A k = QD k Q 1 = Q diag {a k + , a 1 !} Q \ We obtain 



fc=0 



t k (QDQ 

kT 



-i\k 



= Q t ' td ' ag ';'- ff - } Q' 1 = Q diag {«*♦ , e- }Q-\ 

k=0 



(11) 

Thus, to derive the dynamics of the means, we just have to compute the genera- 
tor A, which can easily be done using Table 2. Indeed, we infer from the table that, 
as h — > 0, 

M n (h) = x a x h + 1 x [(1 - ( ai + fa\)h) + HiPiXh] + 2 x (1 - nJ^Xh + o(h) 

= 1 + ((1 - //i)/3iA - ai)h + o(h); 
M 12 (h) = x a x h + 1 x HiPiXh + o(/i) = HiPiXh + o(/i) 

and similarly for M 2 j(/i). From here and (9) we immediately obtain 



7i #i 

^2 72 



7 fe = (1 - Hk)/3kX - a k , 5 k = fi k f3 k X, k = 1,2. 



Now solving the characteristic equation det(A — a I) = for a we find the eigenvalues 
a± of A given by 



<y± = ^ (7i + 72 ± A) , A = \/(7i - 72) 2 + 45i<J 2 

(cf. similar calculations of the threshold parameter for a somewhat different two-type 
model in Section 8.4 of [49]), with the respective (right) eigenvectors 



u± = 



u± 

1 



71 - 72 ± A 



The eigenvalues <r± are clearly different from each other (in any case, this is guaranteed 
by the fact that a + is the Perron- Frobenius root for the quasi-positive matrix A, 
cf. Section A. 8 in [53]), which ensures that A is diagonalizable and we can take the 
transformation matrix Q from (10) to be given by 



Q 



u + u_ 

1 1 



Q 1 = 



U + — M_ 



1 -u_ 
-1 u + 



Hence (8) and (11) imply that 
1 



M(t) 



u + — u_ 
1 

u + — u. 



U + M_ 
1 1 



u + e 
e 



G+t 



e a+t 







1 -U- 









— 1 U + 








&-t g CT +*) 


-t 


u + e° 


-t _ u _ e «+t 



t > 0. 



(12) 
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4 The eventual composition of the population 



As 0+ > er_, it is clear from (12) that a + is the so-called Malthusian parameter of the 
branching process that determines the long-term behaviour of the process means. As 
we said in Section 1, the most interesting for us is the case of supercritical processes 
for which er + > 0, implying unbounded exponential growth of the population (unless 
it becomes extinct at a pretty early stage). Otherwise, the process {Z(t)} would be 
doomed to die out very soon, so that no epidemic would arise. 

It is clear that, in the supercritical case, the ratio of the time t expected number of 
type 2 particles to that of type 1 particles will be given by 



Ri(t) 



M 12 (t) u + u_{ 



M u (t) -u + e CT +'--u_e 



(7 — t 



= \u. 



-At 



t — > oo, 



if the process starts with a single type 1 particle, and by 



M 22 {t) 
M 21 (t) 



u + e 



O-t 



,<T+t 



e a + t 



oO-t 



U. 



l+(l--i + (l)) e 



-At 



t — y oo, 



provided that it started with a type 2 particle (recall that w_ < and a + — er_ = A > 0). 

Therefore, using (7), we see that, regardless of the initial state (z±, z 2 ), the eventual 
ratio of the mean number of T 2 -infected hosts to that for Tx-infected hosts will be one 
and the same quantity 



R = lim 



t— s-co 



E(Z 2 (t)\Z(0) = ( Zl ,z 2 )) 
E(Z l (t)\Z(0) = (z l ,z 2 )) 



\u. 



7 2 - 7i + A 

25,. 



which is a well-known fact from the theory of multitype branching processes (it follows, 
for instance, from Theorem 1 on p. 185 in [5], see also p. 203, ibid.). The convergence 
rate is clearly exponential: the remainder term decays as e~ At . Note also the obvious 
facts that -Ri(O) = 0, R 2 (0) = oo and that R\(t) (R 2 (t)) is an increasing (decreasing) 
function of t (so that always R±(t) < R 2 (t)). 

Thus the single value R = R(a±, a 2 , j3i, 2 , fj,i, /i 2 , X) completely specifies the even- 
tual balance of the mean numbers of individuals of different types in our supercritical 
process (whatever the initial values). This reflects a much deeper result on the long- 
term behaviour of {Z(t), t > 0} — namely, the fact that, with probability one, the 
scaled vector e~ a+t Z(t) will converge, as t — > oo, to a non-trivial random vector whose 
distribution is concentrated on the ray {rv, r > 0} collinear to the (positive) left eigen- 
vector v of A corresponding to the Perron-Frobenius eigenvalue a + (see e.g. Theorem 2 
on p. 206 in [5] and references therein). This implies that convergence to R holds not 
only for the ratio of the means, but for the random variables Z 2 (t) / Zi(t) as well: if we 
denote by A the event {Z(t) ^ for all t > 0}, then 



inn 



R 



on 



A 
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(up to an event of probability zero). In words, this means that either the branching 
process becomes extinct in finite time or the sizes of the subpopulations of individuals 
of the two types grow unboundedly in such a way that their ratio tends to R. 

Observe also that the above shows how fast the composition of the population will 
change if the encounter rate A switches to another value. Suppose that the initial value 
is A'. As we saw, after some (exponentially short) time, the balance of types in the 
process will establish around the value R(X') = R{a\, a 2 , (3 2 , fa, fa, X'). Now if the 
value of A quickly changes to A", then the population will re-establish balance at a new 
level R{X") — again exponentially fast, with the rate characterized by the new value 
of A (provided, of course, that the process will still be supercritical, i.e. <t+ > for the 
new value of A) . 

As we discussed earlier, the increase in the encounter rate A ought to be beneficial 
for the parasite type with higher lethality as it gains more time to spread in the host 
population. Indeed, we have 

OR = (q 2 - ai)(7 2 - 71 + A) 
d\ ~ 25 2 A 

since a 2 > a,\ by assumption and it is obvious that I72 — 7i| < A. Fig. 2 displays 
the dependence of R on A varying in (0,20), for different levels of the lethality a 2 (left 
pane) and mutation rate fa (on the right) of the type 2 pathogen. 




Figure 2: The plots of R as a function of A G (0, 20). For the fixed common values a\ = 0.5, 
fa = 0.2, fii = 0.3 and (3 2 = 0.6, the left pane displays the plots of R for four different 
lethality levels a 2 = 1, 1.5, 2 and 2.5 for fixed fa = 0.2 (the lower the value of a 2 , the higher 
the curve), whereas the right one shows the plots for different mutation rates fa = 0.2, 0.3, 0.4 
and 0.6 (the higher the mutation rate, the lower the curve), for fixed a 2 = 2. 

This example demonstrates that our model is capable of reproducing the growth 
of the frequency of more lethal parasites in a host population when the density of 
the hosts increases. Observe that the threshold value R = 1 (after which type 2 
pathogen dominates the population) may play no critical role: our model is a crude 
approximation for the initial stage of an epidemic only, so Rj [R + 1) will just give a 
proportion of the carriers of type 2 parasite in the population of the infected hosts at 
the end of that stage. Even a relatively law value of that quantity may be fatal for the 
population. 
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It is important to note that it is the "splitting" of the single "virulence" characteris- 
tic of the parasite into two (lethality and transmissibility) that made such a capability 
possible: if, say, there is no difference in lethality (cki = a 2 ) then, as a simple algebraic 
calculation shows, the value of R does not depend on density A. The last observation 
we could have actually made earlier, as it follows from the model construction. 



K,(t) Ri(l) 




012345 012345 



Figure 3: Convergence of Ri(t) to R as t — > 00. For a common set of parameter values, the 
plots display the behaviour of Ri(t), i = 1,2, for A = 2,6, 10 and 14 (from left to right, top 
to bottom), with the respective values R 0.210, 1.076, 1.500 and 1.699. In all the cases, 
the process is supercritical (<r + > 0). 

Figure 3 illustrates the stated exponentially fast convergence of the ratios Ri(t) to 
a common limit R in four situations that have common parameter values cti = 0.5, 
a 2 = 1.5, /ii = /i 2 = 0.2, (3i = 0.3, (3 2 = 0.6, but different encounter rates. The plot in 
the top left corner displays the graphs of Ri(t) < R 2 (t) in the case when the encounter 
rate A = 2 is small enough to allow type 1 parasite to dominate (R ~ 0.210). On the 
top right plot, we see that, for A = 6, there establishes a rough balance (R ~ 1.076), 
whereas on the plots in the second raw we see type 2 parasites to gain dominance 
pretty fast (which is due to higher values of A = a + — <x_), with the limiting values 
R « 1.500 and 1.699, resp. 

The character of the dependence of R on the transmission probabilities is illustrated 
in Fig. 4. For four different values of the encounter rate (A = 2, 6, 10 and 14), the figure 
shows the plots of R as a function of (/3i, fa), restricted to the regions where the process 
is supercritical (i.e. 0+ > 0). The values of the other parameters are a% = 0.5, a 2 = 1.5, 
A*i = A*2 = 0.2. As one could expect, the value of R strongly depends on (fa, fa) and 
is an increasing function of fa and a decreasing one of fa. 

In all the four cases presented in Fig. 4 the threshold value R = 1 is exceeded only 
when the transmissibility of type 2 pathogen is greater than that for type 1 (j3 2 > flx), 
so it may appear that that inequality is a necessary condition for type 2 to prevail. 
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Figure 4: The plots of the limiting value R as a function of the transmissibility probabilities 
(Pi, $2) i n f° ur cases: A = 2,6, 10 and 14 (from left to right, top to bottom), when all the 
other parameters of the model are common (a% = 0.5, «2 = 1.5, fJ-i = fJ.2 = 0.2). The plots 
are restricted to the regions where the respective processes are supercritical. 

This, however, is not true: it turns out that a higher mutation rate from type 1 to 
type 2 can compensate for some lack of transmissibility. Figure 5 shows the plots of 
R as a function of the mutation probabilities (/ii,/^) £ (0, l) 2 for different encounter 
rates, all other parameters being fixed and common, with the transmission probability 
for type 1 being double that for type 2 = 2/3 2 = 0.4). On the left plot corresponding 
to A = 4, the maximum value of R barely exceeds 0.5, whereas on the right one, due to 
the increase in the encounter rate to A = 7, not only the supercriticality region is much 
bigger, but also the maximum value is R — 3. We see type 2 parasite's domination 
in the region where the mutation rate fii (from type 1) is high enough, while /12 is 
relatively small. 

Finally, we turn to the dependence of R on the lethalities a,-. As one can easily see, 



This is quite natural, as the increase in a pathogene type's lethality does not improve 
its chances to prevail when all other parameters of the model remain unchanged. The 
character of the dependence is illustrated in Fig. 6 showing R as a function of («i, a 2 ) £ 
(0, 7) 2 , for A = 2, 6, 10 and 14. Note how the encounter rate A influences the size of 
the region where the process is supercritical (thus, for A2 it shrinks to a narrow strip 
in the («i, a 2 )-domain, corresponding to small values of a\). 




dR 



dR 

da 1 



> 0. 
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Figure 5: The plots of R as a function of the mutation probabilities for A = 4 (the left pane) 
and A = 7 (on the right), for fixed common values of all the other parameters of the model 
{a\ = 1, «2 = 2, Pi = 0.4 and 02 = 0.2. The plots are restricted to the regions where the 
respective processes are supercritical. Note the unusual orientation of the ^-coordinate axes 
(chosen so as to have a better view of the plots). 

5 On the effects of the change in enclosure size on 
the encounter rate 

As we said, our main motivation was to model the effect that the change in "effective 
density" represented by the parameter A (the rate at which infected hosts encounter sus- 
ceptibles) can have on the ratio of the number of hosts infected, say, with T2 pathogen 
to that for Tx-infected hosts. But how does the value of A relate to physically mea- 
surable parameters of the modelled situation — for instance, the density of fish in an 
aquaculture tank (given the tank size and all other parameters of the model are fixed) 
or the size of the tank (given the number of fish in it is fixed)? How will A change 
if one, for example, "squashes" the same host population to a "world" whose linear 
dimensions are twice as small as for the original one? 

The answer to this type of questions will in general depend on what one assumes 
about the character of the hosts' movements (and of course, on the the pathogen 
transmission mechanism — but we will not address this aspect in our simple analysis 
in the section). One of the most popular models for "wandering particles" is the 
famous Brownian motion process {W(t),t > 0} (see e.g. p. 169 in [36]), which can 
be thought of as a continuous analog of a simple (symmetric) random walk. Recall 
that the Brownian motion is denned as a continuous time process with continuous 
trajectories that starts at zero at time t = and has independent Gaussian increments: 
W(t + h) — W{t) ~ iV(0, h) for t,h> 0. One of the key properties of the process is its 
self- similarity: for any a > 0, one has 

{aW(t),t > 0} = {W(aH),t > 0}, (13) 

i.e. these two processes have the same distribution. 

Since the total number of hosts is assumed to be large enough, encounter rates 
are mostly determined by the "local" characteristics (the density of the population 
and the dimensionality of the space) of the enclosure and will have little dependence 
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Figure 6: The plots of R as a function of the lethalities (0:1,02), in the same four cases as 
on Figs. 3, 4: A = 2,6,10 and 14 (from left to right, top to bottom), when all the other 
parameters of the model are common {f}\ = 0.3, = 0.6, ^1 = ^2 = 0.2). The plots are 
restricted to the regions where the respective processes are supercritical and ol\ < ct2 (as 
we assumed). Note the unusual orientation of the a-coordinate axes (chosen so as to have a 
better view of the plots). 

on the "shape of the world" . Therefore, for analysis purposes, we will assume in this 
section that hosts perform independent Brownian motions in a "simple world" S in one, 
two or three dimensions (starting at some "individual" initial points) and show how 
the encounter rates change when one "contracts" the original world without changing 
its shape, i.e. switches from S to the set eS = {ex : x e S}. In this context, the 
dimensionality can actually be though of as a crude description of the the shape of the 
world. 

Suppose that there are N susceptibles in the population and that the movements of 
all the hosts are independent of each other. The value of the parameter A = \n gives 
the average number of encounters of a given infected host with susceptible ones per time 
unit (we assume that n is large enough so that the "conversion" of some susceptibles 
into infected hosts during our modelling "time horizon" does not significantly change 
the number N and hence the encounter rate A). It is clear that A^ = NXi and, 
moreover, that Ai = 1/t*, where t* is the mean time to encounter of our infected host 
with a given susceptible host. Thus the answer to the question on how the encounter 
rate will increase if the host density living in a "fixed world" increases is simple: it 
is just proportional to the number of hosts in the world. However, to understand the 
effect of the world size change (when we switch from S to eS) on the encounter rate 
A = A(e), we will have to analyze that effect on the mean time t* = t*(e). 

The one-dimensional case. First we assume that S = [0, 1] C K. Our hosts move 
according to independent Brownian motions, reflecting from the boundaries and 1 of 
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the set S. This can be formalised by introducing the function 

/ \ f x — I x I if I x I is even, , , , , „ , ^ 

K ' y I — x + \_x\ if [xj is odd, L J L J 

and letting the location of our infected host to be given by Hq = Ho(t) = <p(Ho(0) + 
W (t)) and that of a given susceptible one by Hi = Hi(t) = <p(Hi(0) + Wi(t)), where 
Hi(0) G S are some fixed initial positions and Wi are independent standard Brownian 
motions, i — 0,1. We say that the hosts have an encounter at time t if H (t) = Hi(t). 

Now consider the space eS, where our hosts move according to {sHi(t), t > 0} and 
denote by t*(e) the mean time to encounter of the hosts in this new world, e > 0. It is 
obvious from the self-similarity property (13) that t*(e) = e 2 t*(l) and therefore 

A(e) = e- 2 X(l). 

Thus, if our hosts are confined to a one-dimensional world where they wander at 
random, the encounter rate displays inverse quadratic dependence on the world size: 
say, halving the "living space" will increase the encounter rate fourfold. 

The two-dimensional case. To avoid dealing with any boundaries, we assume that 
our hosts live on the two-dimensional sphere 

§ 2 = {(x, y, z) G M 3 : x 2 + y 2 + z 2 = 1} 

and that our hosts wander on it at random according to independent spherical Brownian 
motions {Hi(t), t > 0} (see e.g. [36], Chapter 15, Section 131), starting at some fixed 
distinct points Hi(0), i — 0,1. 

In this case, Pr(H (t) ^ Hi(t), t > 0) = I, so we need to modify our definition of 
encounter. Fix a small enough 5 > and define the encounter time of the hosts Hi as 
inf{t > : r(H (t), Hi(t)) = 5}, where r(-, •) is the geodesic distance on § 2 . Denote 
by t} = tg(l) the mean value of this time (suppressing the dependence of the initial 
locations Hi(0); to simplify the exposition, we deliberately make it somewhat sketchy). 

Next we consider the "contracted world" eE> 2 , where our hosts wander according to 
{eHi(t), t > 0}, but the definition of encounter remains unchanged (the hosts should 
find themselves within distance 8 of each other); the mean time to encounter for this 
case is denoted by t}(e). Again using self-similarity, we can easily conclude that 

n(e)=e 2 n /£ (l). (14) 

However, we want to relate tg(e) to t|(l), so it remains to clarify the relationship 
between t* s/e (l) and t|(l). 

It is obvious that t* — t*(l) is a decreasing function of rj > 0. As we are interested in 
situations where 5/e is small (despite the small size of the "world", encounters are still 
relatively rare), it suffices to find the asymptotic behaviour of t* as i] — > 0. To do that, 
we first observe that analyzing the dynamics of the position of -f^o(^) relative to Hi(t) 
shows that finding the mean time when the two points are first within distance r\ of each 



19 



other is equivalent to finding the mean time a Brownian particle H*(t) (with an initial 
position at a distance r(H o (0), #i(0)) from the "North Pole" P = (0,0, f) G § 2 and 
local diffusion coefficient V2 times that for the original spherical Brownian processes) 
will need to get within distance r\ of P. 

Denoting by V(t) the projection of H*{t) on the z-axis, one can easily see from Ito's 
formula that {V{t), t > 0} is a diffusion process with the state space [—1, 1] governed 
by the stochastic differential equation 

dV(t) = fi(V(t))dt + a(t)dW (t), 

where W is a standard univariate Brownian motion process and the drift and diffusion 
coefficients are given by 

M*0 = ~2> a \ z ) = ' 

respectively (see e.g. p. 194 in [46]; for convenience we assumed that H*(t) follows a 
standard spherical Brownian motion on S 2 which will have no adverse implications for 
the validity of our analysis). The geodesic distance from H* to P is equal to r\ iff its 
projection on the z-axis equals r := cost], so that we need to find 

v(z) = E(r| V(0) = z), where r = inf{t > : V(t) = r}. 

This can easily be done using the standard technique of the method of differential 
equations (see e.g. Problem B on p. 192 in [36]): the function v(z) is the bounded 
solution to the equation 

1 = -n(z)v'(z) - \o\z)v'\z) ee l -zv'{z) + \{z> - l)v" {z\ z e (-l,r), (15) 

with the boundary condition v(r) = 0. Setting u(z) := (z 2 — l)/4, we notice that 
u'(z) = z/2 and so (15) is equivalent to 

1 = lit) + uv = [uv ) , 

which means that u(z)v'(z) — z + c±, and so 

z + ci A(z + d) 
v (z) = — -— = — . 

V ; u(z) z 2 - 1 

Therefore the general solution to (15) is given by 

v(z) = 2 [(1 + ci) ln(l - z) + (1 - ci) ln(l + z)] + c 2 , 

which is bounded on (— l,r) iff C\ = 1. Now using the boundary condition at z = r to 
find C2 leads to 

1 — 2 

f(z)=41n , zG [-l,r]. 
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To find the asymptotic behaviour of v(z) as rj — > 0, we use cost] = 1 — i] 2 /2(l + o(l)) 
to obtain that, for a fixed initial value z, 

v(z) = 8| I1177I + O(l). 

This means that, for fixed initial positions of H and Hi, we have t* = (c+o(l))\ \nr]\ 
as 7] — > 0. Together with (14) this yields, for small S/e, 
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That is, the encounter rate behaves as 

A(e) « £- 2 M — r A(l). 







\hx6\ - 


\ne\ 



In the two-dimensional case, it might be more natural to relate the encounter rate not 
to the linear dimensions of the enclosure, but rather to its area (proportional to e 2 ). 
The above formula shows that, in this case, the encounter rate in a "shrinking" world 
still grows somewhat faster than the inverse proportional to its area, the latter giving 
the density of hosts. 

The three-dimensional case. We still have (14), and an analysis similar to the one 

carried out in the two-dimensional case shows that now t* = (c + o(l))i]^ 1 as rj — > 
(cf. p.195 in [46]), so that t* s (e) = e 3 t* s (l)(l + o(l)). That is, 

A(,)=^ 3 A(l)(l + o(l)). 

We see that, in the three-dimensional case (assuming that the hosts wander according 
to three-dimensional Brownian motions), the encounter rate in a "shrinking" world is 
inverse proportional to its volume. That is, in this case the rate is proportional to the 
density of hosts. 

To summarise the above analysis, we observe that the encounter rate A = X(e) 
grows rather fast when the linear dimensions (specified by the parameter e) of the 
hosts' "world" diminish. In the three-dimensional case, A is inversely proportional 
to the volume per host, in the two dimensional case it grows slightly faster than the 
reciprocal of the area per host, while in the one-dimensional case the growth rate of 
A is inversely proportional to the square of the size of a host's share of the enclosure. 
This indicates that not only effective density per se, but also the shape of the enclosure 
can be an important factor leading to an epidemic. Thus the nature of the enclosures 
in which animals are kept can be an important factor in determining the progress and 
nature of an epidemic. 



6 A multistage modification of the model 

In this section we will consider an aggregate model for situations in which there are 
several populations of hosts that exist in originally isolated enclosures. 
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There may be many pens with animals at a farm or many fish tanks at an aqua- 
culture facility. Initially, one of the enclosures becomes infected with a single type 
pathogen. This can give rise to a "local epidemic" in the infected enclosure, which 
can be modelled using our processes from Section 2 (assuming that we have super- 
criticality: er+ > 0). The original pathogen may also mutate to become more or less 
lethal. We will initially assume that it may mutate to a more lethal type 2 pathogen 
(a 2 — r oi for some r > 1, a± — a). That new pathogen type can also have different 
transmissibility and mutation rate, but, to make our model as simple as possible, we 
will assume for the time being that it differs from type 1 pathogen in lethality only, all 
other parameters being common. Denoting them simply by (3 and //, we see that the 
Malthusian parameter for that process is given by 

<7+(ai, a 2 ) = X - (2(1 - (j,)p\ - a x - a 2 + a/(«i - a 2 ) 2 + 4/i 2 /3 2 A 2 ) . (16) 

After the epidemic has gone through the initial stage (and so the ratio of the num- 
bers of hosts infected with different pathogen type can be assumed equal to p , where 
p k = R(r k a,r k+1 a, (3, f3, p, p, \)), the infection is transmitted to the next enclosure 
(say, by a worker in a farm situation). The transmitted pathogen is chosen at random, 
so that the probability of transmitting the one with lethality a (denote this event by A) 
is l/(l + po), while the one with lethality ra is transmitted with probability po/(l+Po)- 
This, in turn, may lead to a local epidemic in the new enclosure: we again assume the 
possibility of mutation to a more lethal pathogen (so that now we will have a\ = a, 
a 2 = ra if the event A occurred, and a\ = ra, a 2 = r 2 a otherwise), and to have an 
epidemic we again need + = a + (ai,a 2 ) > (now for the new set of parameters a±, 
a 2 ). Once the epidemic has established itself in the second enclosure (and the balance 
of pathogen types has stabilized around the respective R- value), the next step is the 
transmission of the disease (by means of a random mechanism of the same type as in 
the first instance) to the next enclosure, and so on. 

Scenarios of this type have been encountered often where once a disease is recognized 
in a herd, animals in the infected enclosure are removed or killed, but the disease is 
subsequently found in other herds, for example the spread of Foot and mouth disease 
among herds in Taiwan, which was related to herd size and the number of herds in a 
province [25] . Of course, biosecurity measures are intended to prevent such transmission 
between enclosures, but often the need for diligence is learned after the event. 

It is easily seen from (16) that 

d 

<7 + (0, 0) > 0, — — cr + (a, ra) < 0, lim a + (a, ra) = — 00, (17) 

Oa a^roo 

d 

and also that —R(a, ra) < 0. Thus, if the lethality of the pathogen will keep increas- 
ing, the Malthusian parameter of the model will eventually drop below zero, and then 
the epidemic will collapse. More specifically, setting 

k* = inf{k > : a + {r k a,r k+1 a) < 0}, 
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we see that 

a + (r k a,r k+1 a) < for all k > k*. (18) 

It is clear that the transition of the disease from enclosure to enclosure according 
to the above scheme is described by a discrete time Markov chain {X ra }, the "time" n 
having the meaning of the number of steps (i.e. enclosures infected), X n representing 
the level of lethality of the pathogens in the nth infected enclosure: we set X n = k if 
(011,012) = (r k a,r k+1 a) in the enclosure. Thus the state space of the Markov chain is 
{0, 1, 2, . . . } and the only nonzero entries in the transition matrix P = [Pj,k]j,k>o of the 
chain are 

Pk,k = Pr(-^n+i = k \ X n+1 = k) = — , 

J- + Pk 

Pk,k+i = Pr(X n+1 = k + 1| X n+1 = k) = Pk , 

1 + Pk 

k = 0,1,2,... 

Further, in view of (18), one can assume that once the Markov chain {X n } has 
reached the state k*, the epidemic becomes unsustainable, and hence there will be no 
further transmission of the disease to other enclosures. So we can truncate our state 
space to {0, 1, 2, ... , k*}, which results in a finite decomposable Markov chain with a 
single absorbing state k*. Whatever the current state of the chain, at the next step 
it can either stay at it or move to the right, the transition probabilities forming the 
matrix 

T r T 



M = Q 







where T is the k* x k* substochastic matrix formed by the first k* rows and k* columns 
of P, r — (0, . . . , 0,pk*-i,k*) £ and T denotes transposition. 

The (random) number of steps T the chain will need to reach the absorbing state 
is nothing else but the total number of enclosures that will be affected by the epidemic 
prior to its collapse. Using our model, we can easily find the distribution of T. 

Indeed, using the standard approach to solving such problems (see e.g. p. 80 in [39]), 
we note that as the state k* is absorbing, we have Pr(T < n\ X = 0) = q^l*, where 
are the n-step transition probabilities: 



rpn T 

Vt$\ = Q n = 



n 

1 



rl = (I + T + --- + T n -> T , 



so that for the probability mass vector function f(n) = {fj(n), j = 0, . . . , k* — 1}, 
fj(n) = Pr(T = ra I X = j), one obtains 

f(n) = r n - r n _x = r(T n - 1 ) T , n > 1. 

Of course, we are only interested in the first entry of the vector f(n). 

To compute the mean and higher moments of T one can use the generating function 



/•(*) = 5>"/W = rz{l-zT T y\ \z\ < 1. 

71=1 
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In particular, since j^f*(z) = r(l - zT T ) 1 + rzT T (i - zT T ) 2 , we find that 

= r{l-T T y 2 = l{l-T T y\ 

z=l 

where the last equality follows from the obvious observation that r(J — T T ) 1 = 
= 1 = (1,...,1) G Rf. 
A possible objection to the above simple aggregate model is that pathogens will 
not always mutate to become more lethal. The model can be further generalized by 
allowing, within each enclosure, mutations of our pathogen not only in the direction of 
higher lethality, but also in the opposite direction. So we will first have to generalize 
our basic model from Section 2 to a three-type branching process, assuming that, 
if an enclosure is infected with a pathogen with lethality a, then the pathogen can 
mutate to ones with lethalities r _1 a and ra, respectively, where, as before, r > 1 is 
a fixed number (all other parameters being assumed equal for the different types of 
pathogens). Mathematically, analyzing such processes is essentially equivalent to what 
we did in Sections 2-4, although all the closed-from expressions will become much more 
complicated, and so we will not give much technical detail for brevity's sake. The main 
assertions concerning the asymptotic behaviour of the branching process will remain 
true: there will exist a limiting balance of types in the supercritical case (denote the 
shares of the different pathogens by nj = nj(a), j = 1, 2, 3, J2j 71 j = 1)> which can be 
found from the generator A of the semigroup of the mean matrices, and the almost sure 
convergence of the process scaled by e~ a+t (as before, a + denotes the Perron- Frobenius 
root of A) to a limiting random vector will hold. 

In the multi-stage model, we start with initial infection of one of the enclosures 
with a pathogen with lethality a. That leads to an epidemic (provided, of course, that 
<7+ > 0) in which pathogens of three types will be present, with lethalities given by 
the vector (cci, 0:3) = (r~ 1 a, a,ra). The next enclosure to be infected will receive a 
pathogen chosen at random from those present in the first infected enclosure, and so 
it will have lethality aj with probability ttj, j = 1, 2, 3, and so on. 

Observe that the triplets of lethalities (cti, a 2 , a 3 ) characterizing the pathogens 
present in a given enclosure in our system will all be of the form (x, rx, r 2 x) for some 
x > 0, i.e. lying on a common ray L with the direction vector (l,r, r 2 ). Therefore we 
will again have a basically "univariate" Markov chain {X n } showing what pathogens 
can be present in different enclosures, X n = k meaning that the nth affected enclosure 
was initially infected with the pathogen of lethality r k a, k e Z (and in this case there 
can also be pathogens with lethalities r k ~ 1 a and r k+1 a in that enclosure), assuming 
that Xo = (as a is the lethality of the pathogen that was initially introduced into our 
system). The original state space for the process will be Z = {. . . , —1, 0, 1, . . . }, which 
is infinite in both directions. At each step, the value of the chain can remain unchanged 
(the interpretation being that the pathogen transmitted to the next enclosure had 
the same lethality as the one with which the epidemic started in the current one) or 
can either decrease or increase by one (that is, the transmitted pathogen would have 
lethality values equal to r^ 1 or r times the current one, respectively). 



E(T\X = j) = -t(z) 
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Further, , as before, one can show that J^o~ + (x : rx, r 2 x) < 0, so that, moving along 
the ray L in the "positive direction" , we will eventually enter the subcriticality region 
for the branching process, where the basic reproductive number will be less than one. 
Therefore, at this stage the epidemic will collapse, and hence we again can "truncate" 
the state space for {X n } — now to {. . . , —1, 0, 1, ... , k*}, where k* is an absorbing state 
that has the same meaning as above. 

The variety of behaviours that such a chain can display will be somewhat wider 
than for our first multi-stage model. The dynamics of the chain are determined by the 
behaviour of the mean step values 

E(X n+1 - X n \ X n = k) = n 3 (r k a) - n^a), k < k*. 

In particular, the parameters of the model can be such that the above quantities will 
be negative. Then absorbtion at k* occurs with probability less than one, while on the 
complement event the chain will drift away in the negative direction, which corresponds 
to the disease "fading" , when the pathogen's lethality vanishes, and so on. 

7 Discussion 

We will conclude with a few remarks concerning possible biological interpretation of 
our results. 

The mathematical models we presented show that, at the beginning stages of any 
epidemics that arise in situations where animals (or humans) live in enclosures, the 
density of hosts is an important factor in determining whether more or less lethal 
strains of the pathogen will predominate. The ratio of infections by more and less 
lethal pathogens stabilizes very fast, so that even if measures to prevent further spread 
of the disease are put in place as soon as an outbreak is identified (such as elimination 
of the animals in an enclosure), the relative frequency of pathogen types is likely to 
have changed before action is taken. Other key factors that can also contribute to 
the evolution of more lethal strains and their spread in the host population include 
transmissibility and mutation rates. Mutation rates are known to vary greatly between 
pathogens such as the flu virus, and others such as trematode worms (liver flukes, 
shistosomes etc., which reproduce more slowly). 

Our models show that an increase in the density of animals on farms or mariculture 
facilities will rapidly lead to the dominance of more lethal strains of pathogens if these 
can enter the farm and mutations occur to produce more lethal variants. An example of 
this process has recently been described in the mariculture of fish in [50] . The problem 
has been recognized in intensive poultry production [51]. and the identification of more 
virulent trains of Marburg's disease in chickens [54] may also be an example of this 
process, and it seems possible that the advent of very virulent strains of bird flu during 
the Spanish Flu epidemic may be linked to high densities of soldiers in demobilization 
camps and troop ships, although these men may have been very susceptible due to 
their poor condition [45]. 

Thus the density at which animals are kept should be considered as a risk factor 
for the evolution of more lethal diseases. Assuming chaotic character of movement of 
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animals inside the enclosure where they are kept (as modelled by independent Brownian 
motion processes), we discussed how their effective density affects the key parameter of 
our model specifying the hosts' encounter rate and hence eventually determining what 
pathogen type will predominate. 

The ratio of pathogen types will clearly also depend on the transmissibility of 
pathogen strains. We have focused on differences in lethality between pathogen strains, 
because pathogen strains with increased transmissibility would obviously become more 
prevalent. It is interesting to note that our model showed that, even if the transmis- 
sibility of pathogen of one type is lower than that for the other, the former pathogen 
can still prevail provided it is favoured by mutation. Density of the animal hosts may 
itself increase transmissibility, due to stress on the hosts caused by crowding [52]. In- 
creased transmissibility might in turn lead to the evolution of increased lethality [20]. 
An interesting question is whether increased lethality would be linked in most cases to 
increased transmissibility. Some previous models [20, 21, 17] have assumed that more 
rapid production of copies of the pathogen inside the host would both increase the 
likelihood of transmission to new hosts and also shorten the life of the host. We did 
not make any assumptions of that kind for our model. 

Modern animal husbandry often involves a large number of separate enclosures, 
each containing a large number of animals at very high densities. Once a disease is 
detected in an enclosure, farmers would usually sacrifice or remove the animals, but 
pathogens may be carried between enclosures by various mechanisms, depending on 
the type of pathogen and the biosecurity practices followed). A multistage version of 
our model for this situation suggests that if pathogens are transferred a number of 
times, then the evolution of more lethal pathogens may be very rapid, but the increase 
in lethality will eventually lead to the epidemic becoming unsustainable (hosts dying 
too fast to be able to transit the pathogen). 

We suggest that the outcomes predicted by the mathematical models discussed in 
the present paper can carry important messages for animal husbandry, where there are 
strong commercial incentives to increase the densities of animals in enclosures to very 
high levels, and often very large numbers of enclosures are built in a single farm. 
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